package halberd.utils;

import halberd.data.Complex;

public class FFTUtils {
	public static void fft(Complex[] src, int row, int width, Complex[] dest) {

		Complex[] temp = new Complex[width];
		for (int k = 0; k < width; k++) {
			temp[k] = src[row * width + k];
		}
		temp = fft(temp);
		for (int k = 0; k < width; k++) {
			dest[row * width + k] = temp[k];
		}
	}

	public static Complex[] fft(Complex[] x) {
		int N = x.length;

		if (N == 1)
			return new Complex[] { x[0] };

		if (N % 2 != 0) {
			throw new RuntimeException("" + N + " is not a power of 2!");
		}

		Complex[] even = new Complex[N / 2];
		for (int k = 0; k < N / 2; k++) {
			even[k] = x[2 * k];
		}
		Complex[] q = fft(even);

		Complex[] odd = even;
		for (int k = 0; k < N / 2; k++) {
			odd[k] = x[2 * k + 1];
		}
		Complex[] r = fft(odd);

		Complex[] y = new Complex[N];
		for (int k = 0; k < N / 2; k++) {
			double kth = -2 * k * Math.PI / N;
			Complex wk = new Complex(Math.cos(kth), Math.sin(kth));
			y[k] = ComplexUtils.plus(q[k], ComplexUtils.times(wk, r[k]));
			y[k + N / 2] = ComplexUtils.minus(q[k], ComplexUtils.times(wk, r[k]));
		}
		return y;
	}

	public static Complex[] ifft(Complex[] x) {
		int N = x.length;
		Complex[] y = new Complex[N];
		for (int i = 0; i < N; i++) {
			y[i] = x[i].conjugate();
		}
		y = fft(y);
		for (int i = 0; i < N; i++) {
			y[i] = y[i].conjugate();
		}
		for (int i = 0; i < N; i++) {
			y[i] = ComplexUtils.times(y[i], (1.0 / N));
		}
		return y;
	}

	public static Complex[] dft(Complex[] x) {
		int n = x.length;
		if (n == 1)
			return x;
		Complex[] result = new Complex[n];
		for (int i = 0; i < n; i++) {
			result[i] = new Complex(0, 0);
			for (int k = 0; k < n; k++) {
				double p = -2 * k * Math.PI / n;
				Complex m = new Complex(Math.cos(p), Math.sin(p));
				result[i] = ComplexUtils.plus(result[i], ComplexUtils.times(x[k], m));
			}
		}
		return result;
	}
}
